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ABSTRACT 

TeV-flaring activity with timescales as short as tens of minutes and an orphan TeV flare have been observed from 
the blazar Markarian 421 (Mrk 421). The TeV emission from Mrk 421 is believed to be produced by leptonic 
synchrotron self-Compton (SSC) emission. In this scenario, correlations between the X-ray and the TeV fluxes are 
expected, TeV orphan flares are hardly explained, and the activity (measured as duty cycle) of the source at TeV 
energies is expected to be equal to or less than that observed in X-rays if only SSC is considered. To estimate 
the TeV duty cycle of Mrk 421 and to establish limits on its variability at different timescales, we continuously 
observed Mrk 421 with the Milagro observatory. Mrk 421 was detected by Milagro with a statistical significance of 
7.1 standard deviations between 2005 September 21 and 2008 March 15. The observed spectrum is consistent with 
previous observations by VERITAS. We estimate the duty cycle of Mrk 421 for energies above 1 TeV for different 
hypotheses of the baseline flux and for different flare selections and we compared our results with the X-ray duty 
cycle estimated by Resconi et al. The robustness of the results is discussed. 

Key words: BL Lacertae objects: individual (Markarian 421) - gamma rays: general 
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1. INTRODUCTION 

Mrk 421 is one of the closest (redshift z = 0.03; de 
Vaucouleurs et al. 1991) and brightest blazars known. Due to 
its low-energy synchrotron peak with E sync > 0-1 keV (see, 
e.g., Fossati et al. 2008), it is classified as a high-frequency 
peaked BL Lacertae according to the blazar sequence (Padovani 
& Giommi 1995). Multiwavelength campaigns, especially in 
X-rays (Cui 2004) and y - rays (Tluczykont et al. 2010) have 
shown that Mrk 421 had major outbursts. 21 Moreover, there 
is evidence of correlation between simultaneously measured 

20 Current address: Harvard-Smithsonian Center for Astrophysics, Cambridge, 
MA 02138, USA. 

21 A major outburst usually lasts several months and is accompanied by many 
rapid flares with timescales from tens of minutes to several days, with fluxes 
varying from a few tenths of Crab up to about 10 Crab (see, e.g., Tluczykont 
etal. 2010). 


fluxes in the X-ray and TeV energy band (Fossati et al. 2008), as 
expected within the synchrotron self-Compton (SSC) scenario. 
However, X-ray and very high energy (VHE) emission from 
Mrk 421 do not always correlate (Rebillot et al. 2006) as was 
the case of the TeV flare observed without the X-ray counterpart 
(called “orphan flare”) by Blazejowski et al. (2005). Given the 
limited duty cycle of I ACT instruments, one cannot rule out 
the possibility of lagging counterparts at the other wavelengths. 
Some authors (see, e.g., Reimer et al. 2005 and Sahu et al. 2013) 
have claimed “orphan flares” as evidence of hadronic processes 
taking place in blazars, although non-standard leptonic models 
(see, e.g., Kusunose & Takahara 2006) can also explain them. 

From 2006 to 2008, Mrk 421 was observed by a few 
instruments. For instance, VERITAS and Whipple observations 
from 2006 January and 2008 June do not show significant 
correlations between the y - ray and the optical/radio emission. 
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Moreover, interestingly, a y - ray flare lasting two days was 
detected without increased X-ray activity; unfortunately, the 
data in these wavelengths were not exactly contemporaneous to 
allow the firm conclusion of an orphan TeV flare. MAGIC also 
reported a flare with rapid flux variability in the time period 
2006 April 22-30 (Aleksic et al. 2010). They also detected 
a very intense outburst between 2007 December and 2008 
June that was studied together with simultaneous data in other 
wavelengths. They found that it is difficult to describe the 
spectral energy distribution with the typical variability scale 
of Mrk 421 within the one zone SSC framework (Aleksic et al. 
2012). ARGO-YBJ observed the flux of Mrk 421 to be correlated 
with X-ray emission from 2007 November to 2010 February 
(Aielli et al. 2010; Bartoli et al. 2011). It was pointed out that 
both the X-ray and y-ray spectra harden as the flux increases, 
favoring the SSC model. IACT studies highlight features of 
specific short activity periods of the source, mainly guided by 
external or self-trigger on high states, that could or could not 
be attributed to a general behavior of the source. While the 
sensitivity of Milagro to short-duration flares is less than that 
of IACTs, it is better suited to study long-term variabilities and 
duty cycle, as it operated almost continuously. Mrk 421 was one 
of the brightest sources observed by Milagro and was monitored 
every day for ~6 hr. 

In this paper, we present the analysis of 3 yr of Milagro 
observations (from 2005 September up to 2008 March) of 
Mrk 421 . We provide upper limits on the flux of a flare, limits on 
the flux for different timescales, and an estimation of the y-ray 
duty cycle for energies above of 1 TeV of Mrk 421. 

2. MILAGRO OBSERVATIONS: SIGNIFICANCE MAP 
AND SPECTRUM OF MRK 421 

The Milagro experiment (Atkins et al. 2004) was a large 
water-Cherenkov detector located at 106?68W longitude, 
35?88N latitude in northern New Mexico, USA at an altitude 
of 2630 m above sea level that operated from 2000 to 2008. It 
was designed to detect VHE gamma rays: it was sensitive to ex- 
tensive air showers (EASs) resulting from primary gamma rays 
at energies between 100 GeV and 100 TeV (Abdo et al. 2008a, 
2008b). It had a ~2 sr field of view and a ^90% duty cycle that 
allowed continuous monitoring of the entire overhead sky. The 
main detector consisted of a central 80 m x 60 m x 8 m water 
reservoir with 723 photomultiplier tubes (PMTs) arranged in 
two layers. The top (air-shower) layer (under 1.4 m of water) 
was equipped with 450 PMTs and the bottom (muon) layer (un- 
der 6 m of water) with 273 PMTs. The air- shower layer was 
used to reconstruct the direction of the air shower by measuring 
the relative arrival times of the shower particles across the array. 
The muon layer was used to discriminate between gamma-ray- 
induced and hadron-induced air showers. In 2004, a sparse array 
of 175 “outriggers” was added around the central reservoir. The 
outrigger array covered an area of 40,000 m 2 and each outrigger 
was instrumented with a single PMT. This array increased the 
area of the detector and improved the gamma/hadron separation. 
The instrument reached its final configuration (physical configu- 
ration, analysis procedures, and calibration) in 2005 September. 
This paper only uses data from this last period. 

A detailed description of the Milagro analysis is given in Abdo 
et al. (2012). Here we summarize the information relevant to this 
study. 

Reconstructed Milagro events (hereafter called events) con- 
tain information about the direction (hour angle and declination) 
of air-shower events. From the reconstructed data, sky maps are 


formed. Sky maps are binned in 0? 1 x 0? 1 pixels and contain 
a signal map with the measured counts on the sky and a back- 
ground map with the background expectation calculated using 
the direct integration method described in Abdo et al. (2012). 
The sky maps are constructed for nine independent bins of the 
parameter T (0.2 ^ T ^ 2, in steps of 0.2). This parameter is 
used to give an estimate of the energy of the primary particle 
initiating the EAS and it is defined as 


Nas Nor 

yvUive yylive ’ 
iV AS iV OR 


( 1 ) 


where Nas(or) / N^ e (0R) is the ratio between the number of PMTs 
in the air-shower layer (AS)/outriggers (OR) detecting the event 
and the number of functional PMTs in the air- shower layer 
(AS)/outriggers (OR) at that time. More energetic showers 
contain more particles, cover a larger area, and so fire more 
PMTs: a higher value of the parameter T is then obtained. The 
dependence of the parameter T with the energy of the primary 
particle is shown in Figure 5 of Abdo et al. (2012). 

To maximize the statistical significance when searching for 
sources, a weighted analysis technique is used (Abdo 2007; 
Abdo et al. 2012). A weight is applied to all events, in both 
signal and background maps: gamma-like events are given 
higher weights than cosmic-ray-like events. The values of these 
weights also depend on the parameter T to account for the 
angular resolution of the detector, which is a function of the 
size of the event and of the parameter T. The angular resolution 
ranges from 1?2 for small values of T to 0?35 for large values 
of T (Abdo et al. 2012). To calculate the weights for each 
T bin, the standard Milagro optimization hypothesis for the 
spectrum of an extragalactic source, consisting of a power-law 
with exponential cut-off at an energy of 5 TeV and photon 
index of 2.0, is assumed. This spectrum roughly includes the 
extragalactic background light absorption. The results are not 
strongly dependent on the exact shape of the assumed spectrum. 
For each T bin, the gamma-ray-like excess with respect to 
the background is calculated as the difference between signal 
weighted events and background weighted events. This results in 
nine excess sky maps, one for each T bin. Finally, all nine excess 
sky maps are added into a final excess sky map. The sky map of 
the statistical significance is obtained by using Equation (4) of 
Abdo et al. (2012). 

Mrk 421 was observed with a significance of 7.1 standard 
deviations for a period of 906 days (828 integrated days after 
data quality cuts) from 2005 September 21 to 2008 March 15. 
The median energy of the detected gamma rays is 1.7 TeV 
(under the spectral optimization hypothesis given above). The 
final map of the statistical significance of the excesses in the 
region around Mrk 421 is shown in Figure 1. 

VERITAS has measured the spectrum of Mrk 421 in different 
flux states (classified by level of intensity, from “very low” 
to “very high”; Acciari et al. 2011). In all cases, the energy 
spectrum cuts off below 10 TeV, and an exponential cut-off at 
4 TeV is most typical. A fit to the energy spectrum with Milagro 
data, using the same approach employed to measure the Crab 
spectrum (Abdo et al. 2012), provides a limited constraint to the 
spectrum because the emission from Mrk 421 is concentrated at 
the lowest energy range of Milagro’s sensitivity. Nevertheless, 
we can use the Milagro data to test a specific spectral assumption 
for consistency. As with the Crab measurement, we determine 
the T distribution from the source and generate an expected 
T distribution for several assumed spectra, determining a x 2 
to characterize the agreement between that hypothesis and 
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Figure 1. Sky map of the statistical significance in the region of Mrk 421. The 
significance at the Mrk 421 location (black cross) is 7.1 standard deviations. 

(A color version of this figure is available in the online journal.) 


the data. We find that the VERITAS “low” spectrum is most 
consistent with Milagro data with a x 2 of 12.7 and 9 dof. The 
VERITAS “very low” spectrum is marginally inconsistent with 
the 3 yr integrated average, with a x 2 of 31.1 and 9 dof. The 
“mid” spectrum is inconsistent with a x 2 of 124.1 and 9 dof, 
primarily because of the normalization, rather than the spectral 
shape. Fixing the low-energy spectral index at 2.3 that has been 
measured for Mrk 421 by VERITAS at low-TeV energies, we 
find an exponential cut-off energy between 2.2 TeV and 5.6 TeV 
at 1 standard deviation of confidence, consistent with VERITAS 
measurements. 

3. VARIABILITY 

The light curve (LC) of Mrk 421 is obtained by converting the 
measured weighted event excesses (hereafter called weighted 
excesses) into fluxes through the calculation, with Monte Carlo 
simulations, of the expected weighted excesses for an assumed 
Mrk 421 spectrum. The assumed spectrum is taken with an index 
of 2.3, a cut-off energy of 4 TeV, and a normalization of 0.46 x 


10“ 10 cm -2 s -1 TeV -1 . The spectral index is consistent with the 
VERITAS “low state,” the energy cut-off of 4 TeV is the most 
typical value for Mrk 421 (Krennrich et al. 2001; Aharonian 
et al. 2002, 2005; Konopelko et al. 2008; Acciari et al. 2011), 
and the normalization is obtained by fitting the spectrum as 
described in the previous section. Thus, the measured weighted 
excesses are divided by the expected weighted excesses and then 
multiplied for the integrated flux of the assumed spectrum for 
energies above 1 TeV (0.20 x 10 -10 cm -2 s -1 ). 

The LC obtained with Milagro for energies above 1 TeV is 
shown in Figure 2. Milagro data were recorded on tapes with 
each tape containing data collected over a time interval that, on 
average, is about 1 week. Each time bin in the LC corresponds 
to data recorded in one tape. 

If we assume a constant flux from Mrk_421, we obtain an 
average flux for energies above 1 TeV of Gviiiagro = (0.205 zb 
0.030) xlO -10 cm -2 s _1 . This is consistent with the spectrum 
used to calculate the expected weighted excesses. The x 2 is 134 
for 122 dof, which gives ax 2 probability of 21%, indicating 
that the Mrk 421 flux, measured by Milagro, is consistent with 
being constant during the 3 yr monitoring period. This average 
flux corresponds to (0.85 ± 0.13) Crab, using the Crab flux 
as measured by Milagro (Abdo et al. 2012) for energies above 
1 TeV. 

We also compute the LC of Mrk 421 for energies above 
300 GeV, shown in Figure 3, to make a direct comparison 
with other VHE observations (see Aleksic et al. 2010; Acciari 
et al. 2011). The data from the other instruments have been 
combined 22 to match the Milagro binning. 

As mentioned in the previous section, the spectrum observed 
by Milagro is consistent with the spectrum in the low state 
observed by VERITAS and with being constant over time. 
Therefore, either there are many bright flares that last a much 
shorter time than a week, or there are only a few very bright 
flares such that the average flux over years is still consistent 
with a low state. The fluxes corresponding to the mid and high 
states reported by VERITAS are also shown in Figure 3. IACT 
observations of Mrk 421 during this period indicate that it was 
not in a high state on week timescales and only for one week 


22 The combined average flux has been calculated only by considering the days 
with reported fluxes and assuming the flux to be the same for the whole week. 



Figure 2. Light curve of Mrk 421 (black points) for energies above 1 TeV. The red solid line represents the average value of the flux: ^Milagro = (0.205 ± 
0.030) x 10“ 10 cm -2 s _1 . Each bin represents ~1 week of data. 

(A color version of this figure is available in the online journal.) 
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Figure 3. Light curve of Mrk 421 for energies above 300 GeV. In black: Milagro data; in magenta and in blue: data by Whipple and VERITAS, respectively (Acciari 
et al. 201 1); in green: MAGIC fluxes, calculated as the integral above 300 GeV of the measured spectra (Aleksic et al. 2010). Data from IACT observatories have been 
combined to match the Milagro time-binning. Also shown are the Milagro average flux (1.4 x 10 -10 cm -2 s -1 ) in red, and in gray and cyan the fluxes corresponding 
to the VERITAS mid and lowest high state (“high state A”) of Mrk 421 (3.1 x 10 -10 cm -2 s -1 and 4.1 x 10 -10 cm -2 s -1 , respectively, calculated as the integral above 
300 GeV of the corresponding best fit spectrum with a fixed cut-off energy of 4 TeV; see Acciari et al. 201 1). 

(A color version of this figure is available in the online journal.) 


it was just above the mid state. Clearly, this statement does not 
stand for much shorter timescales than about a week as one can 
notice from the original IACTs LCs (Aleksic et al. 2010; Acciari 
etal. 2011). 

Some of the Milagro measurements correspond to fluxes 
consistent with the high state. However, this is just a result of the 
statistical fluctuations associated with the large error bars of each 
measurement. Thus, it cannot be concluded from Milagro data 
that Mrk 421 was observed in a high state for those bins. In fact, 
all measurements are within 3 standard deviations of the Milagro 
average flux except the two at 53888 and 53958 MJD. These bins 
are above the Milagro average flux at significance of 3.28a and 
3.34a, respectively, but only 1.54a and 1.64a after correcting 
for trials. Therefore there is no significant evidence for flares 
in Milagro data. We can calculate the maximum average flux, 
F max , in a week time period for a flare not to have been detected 
at 99.7% confidence level (C.L.) using the method of Helene 
(Helene 1983). In Figure 4, we show the flux measurements 
given in Figure 3 converted to upper limits (99.7% C.L.) on the 
flux above 300 GeV. The length of the downward arrow for each 
point is the equal to the size of the corresponding error bar of 
each Milagro flux measurement. For comparison, we show the 
Milagro average flux and, the flare observed by VERITAS on 
2008 May with a maximum flux of ~12 x 10“ 10 cm -2 s -1 , as 
reported in Acciari et al. (2011). If we combine 2 the VERITAS 
data to have a weekly time-binning as Milagro, the resulting 
average value is ~5.6 x 10 -10 cm -2 s -1 . 

We have also calculated the largest value of the maximum 
averaged flux (F max ) for flares of different durations. We have 
binned the data in several intervals from one week to six months 
to account for different variability timescales, as outbursts have 
been observed to last up to several months (see, e.g., Tluczykont 
et al. 2010). We have calculated the flux upper limits above 
an energy of 1 TeV. As observed in Figure 5, the values of 
the flux upper limits vary from 2.26 x 10 -10 cm -2 s -1 to 
0.56 x 10“ 10 cm -2 s -1 for a variability timescale of a week 
to six months, respectively. 

HEGRA observed a flare that lasted about three months 
(Aharonian et al. 2003). The maximum flux of the flare ob- 
served in the night of 2001 April 1, was 2.5 x 10 -10 cm -2 s _1 . 


Thus, considering the Milagro upper limits for one (2.26 x 
10 -10 cm -2 s -1 ) and two (1.6 x 10 -10 cm -2 s -1 ) weeks, a flare 
with the maximum flux observed by HEGRA could not have 
lasted longer than one week without having been detected by 
Milagro. 

4. DUTY CYCLE 

As previously mentioned, blazars are highly variable sources 
in short timescales. The lowest steady flux level is called the 
baseline state. The level of activity of a source can be measured 
as the percentage of time that the source spends in flaring states, 
also called duty cycle, given by, 

y\ ti 

duty cycle = DC = ^-L , (2) 

^obs 

where ti is the time that the source spends in the /-flaring state, 
with / running over all the flaring states in the observation period 

^obs- 

To calculate the duty cycle, a flare flux threshold must 
be established to distinguish flaring states. For example, 
Krawczynski et al. (2004) estimated the X-ray duty cycle of sev- 
eral blazars (including Mrk 421) as the fraction of time during 
which the flux exceeds the flux threshold equivalent to 150% of 
the time averaged flux. The same flare flux threshold was used by 
Wagner (2008), with the additional condition that the 50% de- 
viation from the time average flux (considered by Krawczynski 
et al. 2004) was required to be significant at the 3a level. 
Tluczykont et al. (2007) estimated the TeV duty cycle of 
Mrk 421 by using arbitrary flare flux thresholds of 1, 2, 3, 4, 
and 5 Crab. Finally, Resconi et al. (2009) estimated the baseline 
flux R C har and the associated error a C h ar for Mrk 421 in X-rays 
and then calculated the X-ray duty cycle by considering a flare 
flux threshold equal to R Na = R C h ar + Ncr c for N= 1, 3. In 
other words, flaring states are those whose fluxes are 1 and 3 
standard deviation above the baseline flux. 

Choosing a flare flux threshold in terms of the time average 
flux does not allow a direct comparison of duty cycles between 
sources and between different energy bands to be made because 
the time average flux is influenced by the level of activity. In 
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Figure 4. Black arrows: maximum averaged flux in a week time period for a flare in order to have not been detected at a 99.7% C.L.; the length of the downward 
arrow for each point is equal to the size of the corresponding error bar on the flux. Blue solid line: average Milagro flux above 300 GeV. Magenta star: observation by 
VERITAS corresponding to a flaring state; see Acciari et al. (2011); VERITAS data have been combined to have a weekly time-binning. The magenta line marks the 
flux level of the VERITAS flare along the whole Milagro observation period. 

(A color version of this figure is available in the online journal.) 
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Figure 5. Upper limits on the flux as a function of the flare duration. 
(A color version of this figure is available in the online journal.) 


the case of a highly active source, the time average flux, and 
consequently its flare flux threshold, would be much higher than 
the baseline flux, so the duty cycle only refers to the highest flux 
states. On the contrary, for a less active source the time average 
flux is close to the baseline flux, so the duty cycle refers to 
almost all the flaring states. Since the duty cycle of the active 
source includes only the highest flux-flaring states, it is possible 
to obtain a duty cycle value smaller than the one for the less 
active source and erroneously conclude that the latest source 
is more active than the former one. A similar situation happens 
when comparing duty cycle in different energy bands for a given 
source. The case of an arbitrary choice of the flare flux threshold 
will also lead to wrong conclusions, since the same threshold 
selects higher (compared to the baseline flux) flaring states in the 
less active source. Alternatively, the flare flux thresholds defined 
by Resconi et al. (2009) selects flaring states with fluxes equally 
significant (in terms of standard deviations) when compared 
with the baseline flux, independent of the source and of the 
energy band considered. Therefore, we have adopted the flare 
flux threshold definitions proposed by Resconi et al. (2009), in 
order to get an estimate of the TeV duty cycle of Mrk 421 to be 
directly comparable with the X-ray duty cycle. 


The duty cycle definition (Equation (2)) cannot be used 
directly because of the lack of a complete and systematically 
observed set of TeV flux states over a period of years, as in 
the case of X-ray energy band. Instead, we take advantage 
of the time average flux observed by Milagro as follows: 
the total fluence observed by Milagro, Tkmagro x T M ii ag ro (see 
Section 3), where Tkiiiagro is the Milagro observation period, can 
be decomposed into the fluence from the baseline state and the 
fluence from flaring states, i.e., 

^Milagro x ^Milagro = ^baseline x ^baseline + ^ ' ^flar e,z U ? (3) 

i 

where Tb aS eiine is the baseline flux, 7b aS eiine is the time that the 
source spent in the baseline state, Tfl are / is the average flux of 
the /-flaring state in a timescale t i9 with / running over all the 
flares in the Milagro observation period^ 

As discussed in Sections 2 and 3, ^Milagro is constant over 
^Milagro and consistent with the “low” state (not the lowest) ob- 
served by VERITAS. Thus, we could infer that the contribution 
from the flaring states is small. However, this is not necessarily 
correct as it depends on the value of the baseline flux. The con- 
tribution of the flaring states to the total fluence alone does not 
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Figure 6. In red: TeV duty cycle calculated by considering as flares all the states having a flux above 1 TeV greater than Fbaseiine+A^G (left- N = l, right: N = 3); 
the shadowed pink area represents the error associated to the uncertainty on /. In black: X-ray duty cycle, calculated by considering as flares all the states having a 
flux greater that R c har+Ncr c h ar (left: N = 1, right: N = 3; Resconi et al. 2009); the shadowed green area represents its error. The dashed blue line marks the duty cycle 
values for Fb aS eiine = 0.33 Crab. 

(A color version of this figure is available in the online journal.) 


determine the activity of the source, since the same value could 
be obtained by considering many long-duration low-flux flaring 
states or a few short-duration high-flux flaring states. 

-^baseline IS equal tO ^Milagro -^flare? where 7fl are = U fS 

the total time that the source spends in flaring states. If the 
second term in the right side of Equation (3) is rewritten as 
Tflare x (Efl are ), where (Efl are ) is the average flux of flaring 
states, we can solve Equation (3) for and then Equation (2) 
becomes 

(^Milagro ^baseline) 

(Gfl 

are ) ^baseline 

where we have used r M ii ag ro = G obs . 

The estimation of (Ffl are ) is obtained from the distribution of 
day- wise fluxes of Mrk 421 collected by Tluczykont et al. (2010) 
from several VHE experiments (HEGRA, HESS, MAGIC, CAT, 
and Whipple/ VERITAS) from 1992 to 2009. Tluczykont et al. 
(2010) combined the day- wise LCs from different experiments 
by converting the measured flux values to flux values, F , in 
units of the Crab Nebula flux and normalizing to a common 
energy threshold of 1 TeV. This was done by using the energy 
spectrum of the Crab Nebula as measured by each experiment. 
The resulting distribution is well described by a function, g(F), 
which is the sum of: (1) a Gaussian, g G (F), whose mean equal to 
(0.3285 zb 0.0249) Crab (~0.33 Crab) represents the upper limit 
on ^baseline and (2) a log-normal function gi n (F) that describes 
the flaring states (see Figure 3 in Tluczykont et al. 2010): 


8g(F) = 


N g 


og 


\/2n 


exp 


F - flQ 
<* g 


( 5 ) 


and 


gm(F) = 


Ni n 


F cri n 


exp 


(log(-F) - /i] n ) 

2°ln 


2"1 


( 6 ) 


Therefore, we have extrapolated, when needed, the function 
g\r\(F) down to the flare flux threshold F t h r and calculated the 


average flare flux as follows: 


(Glare) = 


f^ m g\ a (F)dF ’ 


( 7 ) 


where F\\ m = 10 Crab is the maximum flux observed in the 
distribution by Tluczykont et al. (2010). (Ffl are ) depends on the 
value of Fthr. For instance, (Ffl are ) is 1.67, 1.84, and 2.64 Crab 
for F thr of 0, 0.33, and 1 Crab, respectively. 

Given Equations (4) and (7), we have calculated the duty 
cycle for several different assumptions of F base iine and for two 
different flare flux thresholds: F thr = Fb ase iine + N x a G , with N = 
1 , 3 (with (Tq defined in Equation (5)). We have also estimated its 
uncertainties as a function of the errors associated with g\ n (F), 
Gum, and Gviilagro • 

The uncertainty in the TeV duty cycle related to the errors 
in the parameters of g\ n (F) has been estimated to be of 4% 
(Patricelli et al. 2013a). The extrapolation of gi n (G) for F\ im 
above 10 Crab is not trivial since it depends on several factors 
such as, e.g., the total available energy of the source and the 
capability to maintain a high flux for a time equal to the 
duration of the flux states considered by Tluczykont et al. (2010). 
Nevertheless, we find that changing Fn m from 10 Crab to 15 
Crab lowers the calculated duty cycle by between 6% and 8% 
depending on the baseline flux (Patricelli et al. 2013b). In the 
following analysis, we do not make further assumptions on F\\ m 
and we only consider the case F lim = 10 Crab. The uncertainty 
on the duty cycle values only considers the error associated with 

Gyiilagro- 

The value of duty cycle, in the case of N = 1 (shown in the 
left panel of Figure 6) ranges from (51 zb 8)% to (32 d= 8)% for 
Gb a seiine = 0 Crab and 0.33 Crab, respectively, while for the case 
of N = 3 (shown in the right panel of Figure 6) it ranges from 
(46 =b 7)% to (27 zb 7)% for F baS eiine = 0 Crab and 0.33 Crab, 
respectively. 

For comparison, the X-ray duty cycle values determined by 
Resconi et al. (2009) are represented by black lines in Figure 6. 
For the case of N = 1, the X-ray duty cycle equal to (40.3 zb 
1.0)% is consistent, within the error bars, with the TeV duty 
cycle almost independent of the value of Fb ase iine- This result 
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could be explained if the X-ray and the TeV activity of the source 
are tightly coupled. However, it should be considered that the 
X-ray duty cycle may be overestimated since fluctuations in 
the X-ray baseline flux have not been discriminated from the 
flaring states, contrary to what we have done for the TeV duty 
cycle calculation by using ginCT 1 ) instead of g(F). The X-ray 
duty cycle, for N = 3, is equal to (18.1 zb 0.5)%, slightly lower 
than the TeV duty cycle independent of the assumed value of 
^baseline; however, the uncertainty in the TeV duty cycle is too 
large to claim a higher activity in TeV than in X-rays. The TeV 
duty cycle becomes consistent with the X-ray duty cycle for 
Tiim >18 Crab. We should note that this result is sensitive to 
other possible emission mechanisms besides the SSC and that 
the TeV duty cycle refers to the 3 yr of Milagro monitoring, 
while the X-ray duty cycle refers to a period of more than 10 yr: 
to do a more direct comparison between the two duty cycle, they 
should be calculated with data collected over the same period 
of time. 

For completeness, we also considered a flare flux threshold 
as given by Krawczynski et al. (2004). In this case, the TeV 
duty cycle cannot be calculated directly from Equation (4) 
since F t h r = F^ = 1.275 Crab corresponds to a flux value 
of standard deviations above the baseline flux (considering 
^baseline = 0.33 Crab). Thus we cannot assume, as done in the 
cases N = 1,3, that the states with flux lower than F th r are only 
fluctuations of the baseline flux. Instead, these states should 
be counted in the total time as well as their flux into the total 
fluence. 

From Equation (2), the ratio between the duty cycle calculated 
with a flare flux threshold F^. and any other flare flux threshold 
F thr is given by 


PCtoE) = Tflare^l) 

DC(Fthr) WfW) ’ 

where 7fl are (F t ^ r ) and T flare (F thr ) are the total time spent in flaring 
states with fluxes above F^ r and F thr , respectively, and are 
proportional to the number of the corresponding flaring states. 
Thus Equation (8) can be rewritten as 


dc(f£) 


DC(F thr ) 

f^g(F)dF 


f 

JkI 


g(F)dF . 


( 9 ) 


Note that the quantity DC(Chr)/ /^ m g(F)dF is independent 

of Fthr, so we can calculate DC(E^) using any of the previous 
estimated values of duty cycle. We then find that the TeV 
duty cycle calculated with a flare flux threshold F^ v ranges 
from (22+ 4 )% to (17^)% for Ebaseiine = 0 and 0.33 Crab, 
respectively. These values must not be directly compared with 
the X-ray duty cycle obtained by Krawczynski et al. (2004) as 
the latter corresponds only to flaring states with fluxes above 
2 standard deviations from the X-ray baseline flux 23 instead of 
9 standard deviations, as considered at TeV energies. The fact 
that the time average TeV flux (as measured by Milagro) is 
much higher than the time average X-ray flux (see footnote 24), 


23 In Krawczynski et al. (2004), the numerical value of the flare flux threshold 
is not reported and cannot be directly estimated, as the average X-ray flux used 
is not given. To get an estimate of this threshold, we considered the X-ray 
average flux estimated by Wagner (2008), equal to 0.86 counts s -1 . With this 
value the flare flux threshold corresponding to the definition of flares of 
Krawczynski et al. (2004) is 1.29 counts s -1 . This value corresponds to 
~^char + 2<r c har, with R c har and <x c har given by Resconi et al. 2009 (0.5 and 
0.4 counts s -1 , respectively). 


when compared with the corresponding baseline fluxes (for 
TeV with the upper limit of 0.33 Crab given by Tluczykont 
et al. 2010, for X-rays see footnote 24), may be an indication 
that the source is more active in TeV energies than in X-rays. 
However, simultaneous observations are needed to make a firm 
statement. 

Finally, using Equation (9), we have obtained a TeV duty 
cycle for a flare flux threshold of 1 Crab that ranges from 27+ 4 % 

to 21+ 5 5 % for Ebaseiine = 0 Crab and 0.33 Crab, respectively. 
These values are lower than (43 ± 13)%, the value calculated 
by Tluczykont et al. (2007), for the ratio of the time in which the 
source was observed in a flaring state and the total observation 
time of the telescopes (IACTs). This is not surprising, as their 
duty cycle may be overestimated because of the observational 
bias of IACTs to continue observations when the source is in 
a high state, leading to an underestimation of the number of 
observations in the baseline state. 

5. CONCLUSIONS 

We have presented results from a 3 yr long-term contin- 
uous monitoring of the BL Lac Mrk 421 with the Milagro 
water-Cherenkov observatory sensitive to gamma rays between 
100 GeV and 100 TeV. Mrk 421 was detected with a statisti- 
cal significance of 7.1 standard deviations over the period from 
2005 September 21 to 2008 March 15. The Milagro measured 
spectrum is consistent with the VERITAS “low state.” Fixing 
the spectral index at 2.3 as measured by VERITAS, we found an 
exponential energy cut-off between 2.2 and 5.6 TeV, also con- 
sistent with the VERITAS measurements. We have also found 
that the Mrk 421 average flux for energies above 1 TeV equals 
(0.205 ± 0.030) x 10 -10 cm -2 s -1 , consistent with being con- 
stant along the Milagro observation period. 

We have found no evidence for flares in the Milagro data. 
Therefore we have established flare flux upper limits for energies 
above 300 GeV for a timescale of ~1 week as a function of time. 
In addition, we have calculated upper limits on the flare flux for 
energies above 1 TeV for timescales from one week up to six 
months, finding that they vary from 2.26 x 10“ 10 cm -2 s -1 to 
0.56 x 10 -10 cm -2 s -1 , respectively. 

Such long-term continuous monitoring has allowed us to 
calculate the y - ray duty cycle of Mrk 421 for flaring states 
with different flare flux thresholds. We have discussed different 
procedures to define the flare flux threshold and justified the 
reasons to adopt the definition given by Resconi et al. (2009) 
in our analysis. Two cases are presented in detail: flare flux 
threshold of 1 and 3 standard deviations above the baseline flux. 
We have compared the corresponding results (see Figure 6) 
with the X-ray duty cycle estimated by Resconi et al. (2009). 
We find that the TeV duty cycle is consistent with the X-ray 
duty cycle and therefore with the SSC emission mechanism, 
although it is sensitive to alternative emission processes. More 
observations and further studies, for instance, of the expected 
correlation between the activity at TeV energies and X-rays, are 
required to reduce the uncertainties in the quantities involved in 
the duty cycle calculation and to obtain a conclusive result on 
the emission mechanisms involved. 

The High Altitude Water Cherenkov (HAWC) detector, the 
successor of Milagro, will be able to produce a more accurate 
analysis of the TeV emission from Mrk 421, with its greater 
sensitivity (10-15 times better than Milagro). In particular, with 
HAWC it will be possible to determine with greater accuracy 
the average flux, as well as the distribution of flux states of 
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Mrk 421, allowing a more precise estimation of the TeV duty 
cycle. 
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